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Relic gravitational waves (RGWs) generated in the early Universe form a stochastic GW back¬ 
ground, which can be directly probed by measuring the timing residuals of millisecond pulsars. In 
this paper, we investigate the constraints on the RGWs and on the inflationary parameters by the 
observations of current and potential future pulsar timing arrays. In particular, we focus on effects 
of various cosmic phase transitions (e.g. e"’"e“ annihilation, QGD transition and SUSY breaking) 
and relativistic free-streaming gases (neutrinos and dark fluids) in the general scenario of the early 
Universe, which have been neglected in the previous works. We find that the phase transitions can 
significantly damp the RGWs in the sensitive frequency range of pulsar timing arrays, and the upper 
limits of tensor-to-scalar ratio r increase by a factor 2 for both current and future observations. 
However, the effects of free-steaming neutrinos and dark fluids are all too small to be detected. 
Meanwhile, we find that, if the effective equation of state w in the early Universe is larger than 1/3, 
i.e. deviating from the standard hot big bang scenario, the detection of RGWs by pulsar timing 
arrays becomes much more promising. 

PACS numbers: 04.30.-w, 04.80.Nn, 98.80.Cq 


I. INTRODUCTION 

Inflation is the most popular scenario of the extremely 
early Universe [1-3] . In addition to elegantly solving the 
flatness puzzles, the horizon puzzles, and the monopole 
puzzles in the hot big bang universe, inflationary models 
predict the primordial density fluctuations (scalar per¬ 
turbations) with the nearly Gaussian distribution, and 
the nearly scale-invariant power spectrum [4] , which have 
been strongly supported by the recent observations on 
the anisotropies of cosmic microwave background (CMB) 
radiation [5-7], and the distributions of the galaxies in 
various large-scale structure observations. 

In the inflation scenario, a stochastic background of 
relic (primordial) gravitational waves (RGWs) in full fre¬ 
quency range was inevitable produced due to the supera- 
diabatic amplification of zero-point quantum fluctuations 
of the gravitational field [8], which provides the unique 
window to study the physics in the early Universe be¬ 
fore the big bang nucleosynthesis (BBN) stage. Nowa¬ 
days, the detection of RGWs in different frequencies has 
been carried out by different methods. In the lowest 
frequency range with / < I0“^^Hz, it can be detected 
by their imprints in the GMB anisotropies in temper¬ 
ature and polarizations [9]. The recent Planck observa¬ 
tions on the GMB temperature and E-mode polarizations 
give the tightest constraint on the tensor-to-scalar ratio 
r < 0.11 [7], which is consistent with the result r < 0.12 
obtained from the joint analysis of BICEP2/Keck Ar¬ 
ray and Planck B-mode polarization data [10]. In the 
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near future, one anticipates that the detection limit can 
arrive at r = 0.01 for various ground-based or balloon- 
borne experiments (including Keck/BIGEP3, POLAR- 
BEAR, SPT-3G, AGTPOL, GLASS, QUBIG, QUIJOTE, 
QUIET, Simons Array, EBEX, PIPER, SPIDER et ah). 
While for the space-based missions, such as GMBPOL, 
LiteBird, COrE, PRISM, PIXIE, the detection limit can 
be r = 0.001 [11]. 

In the high frequency range, i.e. / G (lO”'^, 10^)Hz, 
RGWs are detected by the ground-based and space-based 
laser interferometer gravitational-wave observatories. So 
far the most stringent bound is obtained by the joint 
analysis of the 2009-2010 LIGO and Virgo data, which 
gives energy density of RGWs Ugw(/) < 5.6 x 10“® at fre¬ 
quency band spanning 41.5-1726Hz [12]. In the future, 
the detection limits at / ^ lOOHz are expected to be 
~ 10“® for the second-generation gravitational-wave de¬ 
tectors, such as AdvLIGO, AdvVirgo, KAGRA and so on 
[13]. Eor the third-generation detectors, such as Einstein 
Telescope, this limit could arrive at ~ 10“^^ [14]. The 
space-based detectors are sensitive to the lower-frequency 
gravitational wave (about 0.1 mHz to 1 Hz). Eor the fu¬ 
ture eLISA/NGO mission, the detection ability is hopeful 
to be ~ 4 X 10“^° at / = 4 mHz [15]. In addition, the 
limit can be even improved by 5-10 orders by BBO [16], 
DEGIGO [17] and ASTROD [18] in the far future. 

In this paper, we shall focus on the RGWs in 
the median frequency range. In the frequency band 
/ G (10“®, 10“^)Hz, the gravitational waves are probed 
through the pulsar timing arrays (PTAs). It is well 
known that, the millisecond pulsars are the very stable 
clocks. The tiny timing residuals of these pulsars are 
caused by some intrinsic or environmental noises, as well 
as gravitational waves [19]. It was realized that the tim- 
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ing residuals from an array of pulsars could be analyzed 
coherently to separate GW-induced residuals from other 
effects [20, 21]. So, it provides the excellent way to detect 
RGWs in the median frequency range. Nowadays, these 
observations are carried out by three groups (PPTA, 
EPTA, NANOGrav). Recently, all teams released their 
latest timing results [22-24], and the tightest constraint 
is obtained by the PPTA team, they placed a 95% upper 
limit on the strain amplitude (at a frequency of yr“^) 
in the power-law model of Ag^ < 1.0 x 10“^^ for spec¬ 
tral index —2/3 [22]. These bounds could be significantly 
improved by the potential observations of future Interna¬ 
tional Pulsar Timing Array (IPTA), Five-hundred-meter 
Aperture Spherical Radio Telescope (FAST) in China 
and the planned Square Kilometer Array (SKA) projects 
[25]. These upper limits of GW background have also 
been applied to constrain the RGWs and on the infla¬ 
tion physics [23-28]. In all these previous works, a sim¬ 
ple evolution model is used to calculate the RGW power 
spectrum, where the effects of free-streaming fluids and 
various cosmic phase transitions are neglected. How¬ 
ever, it was known that, the relativistic free-streaming 
gas (such as the neutrinos in the Universe) gives rise to 
an anisotropic term in the evolution equation of gravi¬ 
tational waves, which can significantly damp the RGW 
spectra at / > 10“^®Hz [29, 30]. Another effect is 
caused by the successive changes in the relativistic de¬ 
grees of freedom during the radiation-dominant stage, i.e. 
the QCD transition, e+e“ annihilation, the electroweak 
phase transition and so on. During these transitions, the 
evolution of scale factor was altered compared with the 
standard radiation-dominant stage, which left the im¬ 
prints in the RGW spectrum at / > 10“^°Hz [31-33]. 
So, both effects can change the RGWs at the median 
frequency range. In this paper, we shall investigate in 
details the effects of free-steaming gas and cosmic phase 
transitions on the RGW spectrum in a general cosmo¬ 
logical scenario. In particular, we shall focus on their 
influences on the detection of RGWs by the current and 
future PTAs and the corresponding cosmological impli¬ 
cations. 

The outline of this paper is as follows. In Sec. II, 
we introduce the RGWs in the accelerating Universe by 
considering the effects of free-steaming gas, cosmic phase 
transitions and unusual equation of state. In Sec. Ill and 
Sec. IV, we discuss the constraints on the inflationary 
parameters by the current and future PTA observations. 
Section V summarizes the main results of this paper. 


II. RELIC GRAVITATIONAL WAVES 


The action of gravitational wave hij is [34] 

r _r— 0 ^“^ 1 

S= I dTcPx^f^ , (1) 


where g^iy=diag{—a^, a^, a^, a^} is Friedmann-Lemaitre- 
Robertson-Walker metric, with a the scale factor, G the 


Newtonian gravitational constant and r the conformal 
time. The term H^ is anisotropic stress, which includes 
the contribution of large-scale magnetic field [35], free- 
streaming relativistic particles (e.g. neutrinos after their 
decoupling) [29, 36] and so on. 

By applying Fuler-Lagrange equation to (1), we can 
obtain the equation of motion for /iy (x, r). Decompose 
the perturbation hij (x, r) and H^- (x, r) by Fourier trans¬ 
formation, then we obtain the equation of evolution for 
the mode with conformal wavenumber k: 


h^,{T) -f 2—h'j,(r) +k^hk{T) = 167rGa^(T)nfc(r), (2) 

a 

where the prime (') denotes derivative with respect to t 
and Hfc is the Fourier component of I\.ij{x,T). 

The power spectrum of RGWs at r is defined by [34] 

Ph{k,T) = 647rG^|hfc(T)|^, (3) 

while the initial power spectrum (spectrum at r = 0) is 
usually parameterized by a power-law form [37] 

( ]c \ '^t 

k^) ’ 

where r is tensor-to-scalar ratio, and An^ko) is the value 
of scalar power spectrum at k = /cq, with ko a pivot 
number, and nt the power index of tensorial spectrum. 
Throughout this paper, we use feg = 0.002 Mpc“^, and 
4l/i(fco) = 2.371 X 10“® (this result is obtained by con¬ 
verting An = 2.139 X 10“® [7] at fc = 0.05 Mpc“^ into 
that at fc = 0.002 Mpc“^). 

The tensor-to-scalar ratio r is currently constrained to 
be r < 0.11 (95%.C.L.)[7]. For de-Sitter inflation, the 
power spectrum becomes flat when the modes leave the 
horizon, so n* = 0 [38]. For slow-roll inflation, due to 
the null energy condition, nt < 0 [39], but non-canonical 
inflation or the extension of slow-roll model may predict 
nt>0 [40]. 

It is convenient to convert power spectrum into dimen¬ 
sionless spectrum of energy density [34] 


Dgw(fc, T ) 


1 dpg^ 

Pc dink 


k^Ph{k,T) 

l2a?(T)m(Ty 


( 5 ) 


where pgw is the energy density of gravitational waves, 
and Pc = 3H^/ (87rG) is the critical energy density of the 
Universe (we set c = 1). This expression is equivalent to 
Eq. (7) in our previous work [25]. 

Throughout this paper, we use the standard model of 
cosmology. The cosmological parameters adopted in this 
paper are listed here: Dm = 0.308, Zeq = 3365, cur¬ 
rent conformal time tq = 1.41 x 10'’’ Mpc, present scale 
factor a(ro) = 1, current Hubble parameter H{tq) = 
67.8km s“^ Mpc“^, and the reduced Hubble parameter 
h = i7/(100km s“^Mpc"^) [7]. 

Because the equation of evolution (2) is difficult to be 
solved analytically (see [41], for example), it is useful to 
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introduce transfer function T{k,T) as follows [42] 


T{k,T) 


hkjr) 

/ifc(O)’ 


( 6 ) 


where /ifc(O) is the initial amplitude of the mode with 
wavenumber k. Now all the effects of evolution in the 
post-inflation stage are attributed to the transfer func¬ 
tion and we split the transfer function into several parts 
to account for different damping factors [34, 43] 


T{k,T) = T^(fc, r)xTpT(A:,r)xTj,(/c,r)xTEos(fc,T), ( 7 ) 


where represents the transfer function of the cosmic 
expansion (the reason of cosmological redshift z), while 
Tpx is for cosmic phase transitions in hot Universe, Ti, is 
for free-streaming particles(e.g. neutrinos ly) and TeoS is 
for unusual equation of state (EoS) before BBN epoch. 

Inserting (3), (4) and (6) into (5), the present spectrum 
of energy density becomes 


Ug^(/c, To) 


rk'^Aii{ko) / k Y* 

12H^to) \ko) 


xT\k,To). 


( 8 ) 


From (7) and (8), we could see that the transfer func¬ 
tions determine the strength of RGW spectrum to some 
extent and that they are closely related to the physics 
after inflation. Now we would like to discuss the transfer 
functions one by one. 


A. RGWs in the accelerating Universe 


RGWs were stretched out of the horizon during the 
inflationary epoch and re-entered the horizon in the fol¬ 
lowing epochs of decelerating expansion and recent ac¬ 
celerating expansion [3, 38]. 

When stretched out of the horizon, RGWs started to 
freeze out and the amplitude kept constant, because out 
of the horizon there are no anisotropic tensor to damp 
the RGWs even when there are some particles whose free 
path may be comparable to the horizon [3, 29]. When 
the Universe exited inflation, the horizon started to in¬ 
crease and the modes with high k (thus short wavelength) 
started to re-enter the horizon. The criterion of horizon 
crossing is a(T)H{T) = k [44]. After re-entering the hori¬ 
zon, RGWs are mainly damped by the expansion of the 
Universe. According to (2), when the anisotropic term 
Hij is neglected, the amplitude of modes in the inner part 
of the horizon {a[T)H{T) k) would satisfy 


hk{T) 


fefc(O) 

a(r) ■ 


( 9 ) 


So the evolution of RGWs is closely related to the history 
of expansion. Turner et al [42] formulated the transfer 
function in the decelerating expansion by fitting it to the 
numerically integrated results, while Zhang et al noticed 
the effect of accelerating expansion [45, 48]. Now the 


transfer function in the expanding Universe is the prod¬ 
uct of Trm (radiation and matter dominated era) and Ta 
(A dominated era), or precisely [42, 45-48] 


Tz{k,To) = 


3Umji(feT-o) 

krn 


'1.0-f 1.36-^-b2.50(-^)2, 

^eq ^eq 


( 10 ) 

where keq = 0.073U„iIi^ Mpc ^ is the wavenumber cor¬ 
responding to the mode that entered the horizon at the 
equality of matter and radiation, and ji (kro) is the spher¬ 
ical Bessel function of the first kind. 

In Fig. 1, we plot the energy density spectrum of 
RGWs by only considering the damping effect of the 
cosmic expansion, where different spectral indices rit are 
considered. For the comparison, we have also plotted the 
current and potential constraints/detection of RGWs by 
various observations. 


B. Effects of cosmic phase transitions: e''"e 
annihilation, QCD transition and SUSY breaking 


When the hot Universe gradually cooled down in the 
radiation-dominant era, massive particles became less 
relativistic and contributed less to the energy density 
of radiation. Typical cases are the breaking of super- 
symmetry (SUSY) and the combination of quarks into 
hadrons (phase transition of quantum chromodynamics 
i.e. QCD phase transition). Besides that, the annihila¬ 
tion of electrons and positrons also induced big change 
in the energy density of radiation. Here, we generally 
call all these three effects and similar effects that could 
induce the change of energy density of radiation cosmic 
phase transition. 

The cosmic phase transitions change the energy density 
of radiation thus the Hubble parameter H{t) = a'la^, 
and affect the amplitude of RGWs through Eq. (9). The 
transfer function of phase transitions is [31, 32] 



where 5 *(T'fe) is effective degree of freedom (here effective 
degree of freedom is obtained by converting all particles 
into effective photons, see Chapter 3.3 in [2]) for total 
energy density and g*s(7fe) is the effective degree of free¬ 
dom for entropy density at temperature Tfe, while (?*o 
and (?*so are the present values. The parameter Tk is the 
temperature at which the mode k re-enters the horizon. 

By use of the condition of horizon crossing aH = k 
and the conversation of entropy = constant, we 

can convert temperature into wavenumber: 



( 12 ) 
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FIG. 1: The energy density spectrum of RGWs r2gw(/) for different inflationary models, where we only consider the damping 
effect caused by the cosmic expansion in the standard hot big-bang Uiriverse. In this figure, we also show the current constraints 
on RGWs by different observations (red lines) and the potential coirstraints by future observations (blue lines). The details of 
these constraints can be found in the main text. 


where T(to) = 2.7255 K is the current CMB temperature 
[49] and Zeq is the redshift at r2m('req) = I2r('req). We can 
further convert wavenumber into frequency: k = 27r/. 

The transfer function of phase transitions depends on 
not only the temperature but also the physical charac¬ 
teristics of particles, such as mass, degree of freedom and 
the statistical features (bosons or fermions) and so on. 
We follow the list of particles in [31] but update the data 
according to [50]. The temperature when particles de¬ 
couple from the others is also very important. Here we 
adopt instantaneous decoupling of neutrinos from the ra¬ 
diation at T = 1.5 MeV, and that electrons annihilate 
with positrons immediately at T = 0.1 MeV [51]. The 
phase transition of QCD is treated as first order transi¬ 
tion at 155 MeV [52]. 

Fig. 2 shows the square of transfer function rprp(fc, tq) 
at different temperature. When the cosmic phase tran¬ 
sitions are considered, the energy density spectrum of 
RGWs can be reduced by about 60% at 10^ ~ 10^ MeV 
and by about 20% at 0.1 ^ 10 MeV. When we include 
the SUSY particles, the spectrum could be damped as 
much as 70% at T > 10® MeV. However, SUSY breaking 
is beyond the detecting ability of PTAs (52 ~ 3.9 x 10® 
MeV). So we will not consider SUSY particles later. 



LogioTfMeV] 


FIG. 2: The dependence of rpq.(fc,ro) on temperature: The 
solid line (black) is the case without SUSY particles, while 
the dashed one (red) is the case containing them. On the 
dashed line, from right to left there are three typical steps, 
which are the evidence of breaking of SUSY (~ 10® MeV), 
phase transition of quarks (~ 155 MeV) and the annihilation 
of e'’"e“('~ 0.1 MeV). The figure is plotted without considering 
the dark fluid (introduced in section IIC 2). 
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C. Effects of relativistic free-streaming gases: 
neutrinos and dark fluids 

1. Neutrinos 

Now, let us discuss the anisotropic term in (2). As the 
cosmic magnetic filed is very small, we will ignore the 
magnetic sources [35]. When neutrinos decouple from 


the rest of radiation, the free streaming of neutrinos con¬ 
tributes to the anisotropic stress nfe(r), which plays the 
role of friction in the equation of motion and damps the 
amplitude of RGWs, so it is necessary to study the damp¬ 
ing effect of neutrinos. 

In very high precision, the transfer function of rela¬ 
tivistic free-streaming particles is [34, 36] 


T^{k,To) 


15(324135000 - 48118000/^ -f 3152975/2 - 55770/^ -h 14406/^) 
343(15 + 4/,)(50 + 4/,)(105 + 4/,)(180 + 4/,) 


(13) 


where fi, is the fraction of free-streaming particles’ en¬ 
ergy density over the total energy density. This formula 
applies for modes satisfying ^/2k kgq- When k de¬ 
creases and becomes comparable with fcgq, this formula 
does not apply any more. However, as k decreases, the 
transfer function will gradually increases to 1 [53]. 

Because the sensitive band of PTAs is 10“® ~ 10“^ Hz, 
which is much higher than fceq, we would not consider the 
behavior of free-streaming neutrinos near fceq, thus (13) 
is enough for us and we take it as the transfer function 
of free-streaming particles. 

When there are neutrinos only, is the fraction of 
energy density of free-streaming neutrinos. Before their 
decoupling from the radiation, = 0, thus ry(fc, rg) = 
1, so neutrinos could not damp the RGWs. After the 
decoupling, f^, started to evolve with the Universe and 
depends on the temperature. 

In Fig. 3, we plot the fraction of neutrinos’ energy den¬ 
sity. We find that can be roughly divided into two 
constant stages = 0.491 with T^{k, tq) = 0.589 (before 
e+e" annihilation) and fi, = 0.405 with T^{k, rg) = 0.645 
(after e^e“ annihilation). Therefore, our analysis is con¬ 
sistent with that in [29, 34] 

It should be noted that the dip at logig(T/MeV) £ 
(—1.2,—1.0) and the spike at log^g^T/MeV) £ 
(—1.0,—0.6) are artificial signatures. They are due to 
the imperfect deal with annihilation (we assume instanta¬ 
neous annihilation here, which is not the real situation). 
The same artifacts appears in [31], while [54] explains 
this in its Appendix C. 

Fig. 4 plots the transfer function of cosmological red- 
shift, neutrinos and cosmic phase transitions. Appar¬ 
ently, Tz dominates all the damping effects in a wide 
range from 10“^®Hz to lO'^Hz, while neutrino streaming 
would decorate the spectrum in the low frequency band 
(< 10“^°Hz) and the cosmic phase transitions modifies 
RGWs in the high frequency range (> lO^^^jjz). As we 
could see, phase transitions show their importance espe¬ 
cially around 10“® Hz, which well resides in the sensitive 
band of PTAs. 


2. Dark fluids 

In the Standard Model of particle physics, there are 
three species of neutrinos. But when we consider the 
non-instantaneous decoupling of neutrinos, the effective 
number of neutrino species Nes = 3.046 [55]. The cur¬ 
rent constraint from GMB and matter power spectra is 
A^eff < 3.376 (95% confidence) [56], while the BBN ob¬ 
servations give A^efi < 3.41 (95% confidence) [57], so 
tension is still possible between theory and observation. 
Weinberg [58] proposed massless bosons to relax the ten¬ 
sion and found AA^eff = 0.39 for bosons decouple at 
T > lOOMeV. Thus the current GMB observation can 
not exclude the existence of extra neutrino-like particles. 

Here, we assume the existence of extra fermions to ad¬ 
dress the AAgff problem. Since massive particles freeze 
out at low temperature and do not contribute to NeS, we 
would like to assume the particles are massless. Due to 
the constraints on Aeff, it is well justified to set the physi¬ 
cal degree of freedom of dark fluid to be two (one for par¬ 
ticle, the other for anti-particle), then AA’eff = 1. To con¬ 
sider the possibilities of different particles, we shall con¬ 
sider different decoupling temperature. Since the PTA 
band / £ (10“®, 10“^) Hz corresponds to 52 ^ 3.9 x 10^ 
MeV, we consider three decoupling temperature 10,10^ 
and 10^ MeV in our discussion. As the extra particles can 
be any possible particles beyond current observations, we 
simply call them dark fluids (DFs). 

Fig. 3 shows the fraction of free-streaming particles in¬ 
cluding dark fluid. When we include DF, becomes 
the fraction of all free-streaming particles’ energy den¬ 
sity. Both neutrinos and dark fluid contribute to . Rich 
features occur in the fraction f^, of free-streaming parti¬ 
cles. In Fig. 3, we consider three different decoupling 
temperature (T = 10,10^ and 10^ MeV) of dark fluid. 
Besides the two stages fi, ~ 0.476 with T/(fc, rg) = 0.599 
and fi, ~ 0.563 with r/(fc,rg) = 0.547 separated by 
annihilation of e“'"e“ at T = 0.1 MeV, one more stage 
{fv ~ 0.13, r/(fc, Tg) = 0.867) appears when dark fluid 
decouples at 10 or 10^ MeV, and a third stage {fi, « 0.025 
with T^{k, Tq) = 0.973) occurs when dark fluid decouples 
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FIG. 3: The fraction fi, of free-streamings’ energy density over 
the total energy density: The black curve is the case without 
dark fluid, while the red, green and blue ones contain dark 
fluids which decouple at 10, 10^ and 10® MeV respectively. 
Neutrinos start to decouple at l.SMeV and e''"e“ annihilate at 
0.1 MeV, while phase transition of QCD happens at 155 MeV. 
When there are neutrinos only, the annihilation at T = 0.1 
MeV splits fi, into two stages fi, « 0.491 and 0.405. When 
dark fluids are included, fi, increases and has more stages: 
/„ « 0.13 for T G (1.5,155) MeV and /„ « 0.025 for T G 
(155,1000) MeV. 
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FIG. 4: The compilation of transfer functions from Tz{k,To), 
Tv{k,To) and Tpt(A:,to). The grey line is for Tz, while the 
red one is for the product of Tz and and the dotted blue 
one combines all the three damping factors. The shaded area 
denotes the sensitive band of pulsar timing arrays. No SUSY 
particle is considered here. 

at 10® MeV. Thus dark fluid would damp the RGWs at 
higher frequency. 

By use of the expression of spectrum (8), the trans¬ 
fer function (10), (11) and (13), we can plot the energy 
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FIG. 5: The energy density spectrum of RGW including 
Tzik, to), TpT{k, to) and Tv{k, t). In this figure, r is chosen to 
be 0.1 and nt — 0. The top grey line is the spectrum consid¬ 
ering the effect of cosmological redshift Tz(k,To) only, while 
the pink curve includes free-streaming neutrinos T^ik, tq) and 
cosmic phase transitions TpT{k,To). Besides damping from 
Tz(k,To), TpT{k,To) and Tv{k,To), the black, green and the 
dashed blue line also include the effects of dark fluid, which 
is assumed to decouple at 10,10® and 10®MeV respectively. 
The shaded area denotes the sensitive band 10~® ~ 10“^Hz 
of PTAs. (In this plot, we do not consider SUSY particles.) 


density spectrum of RGWs for given r and nt- 

Fig. 5 shows the spectrum of energy density for r = 0.1, 
nt = 0. In this plot we could see the damping effects 
caused by cosmological redshift, the cosmic phase tran¬ 
sitions and the free-streaming neutrinos/DF. Neutrinos 
damp the spectrum at 10“®® < / < 10“®®Hz, beyond 
the sensitive band of PTAs, thus PTAs cannot detect 
any relic signals from neutrinos. The cosmic phase tran¬ 
sitions mainly damp the spectrum at / > 10“® Hz by 
amount of 60%. The big jump at / ^ 10“® Hz is the 
signature of QCD phase transition (T = 155 MeV). 

In Fig. 5, we find that DFs have three kinds of sig¬ 
nature in the spectrum. Firstly, steps occur after their 
decoupling from the radiation. Secondly, the spectrum 
at / < 10“®° Hz is damped more. Lastly, spectrum at 
/ > 10“® Hz is raised. 

At the moment of decoupling (T = 10,10^, 10® MeV), 
big jumps occur, as the free-streaming dark fluids start 
to damp the spectrum. In the range of the decoupling 
of neutrinos (10“®® ~ 10“®® Hz), the appearance of 
dark fluid increases the fraction of free-streaming par¬ 
ticles thus damps the spectrum more. This result is 
consistent with our intuition. But at higher frequencies 
(/ > 10“®Hz) dark fluids raise the spectrum when com¬ 
pared to the case without dark fluids. 

The reason why dark fluids raise the spectrum at high 
frequencies resides in the massless feature of dark fluid. 
We can clearly see this from the transfer function of phase 
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transition (11). Massless particles can keep relativistic at 
low temperature and always contribute to the effective 
degree of freedom for energy and entropy density. So the 
addition of dark fluids increases the value of g*so and 
g*s{Tk). As g^,{Tk) Ri g*s{Tk) for all Tk (see Chapter 3.3 
in [ 2 ]), then 

TpT(fc,ro)« (14) 

At high temperature (which corresponds to high fre- 
quency) 5*3 (T^) ^ g*sO, so the increase of 5 *s 0 dominates 
and makes Tpx increase a bit. 

D. Effects of general equation of state 

After inflation, the Universe should experience a period 
of reheating to create the matter and dark matter we see 
today. Some time later, dark matter should decouple 
from the rest of particles [3]. It is suspected that this 
epoch before BBN can be dominated by massive particles 
and cause deviation from ideal fluid [34, 59], thus unusual 
EoS may occur. EoS could deeply change the shape of 
spectrum of energy density through scale factor as we 
would see below. 

The relation between scale factor and conformal time 
could be generally parameterized by 

a oc (15) 


transfer function of cosmic expansion Tz(k,TQ) and that 
of unusual EoS 7Eos(fc,7o) into the total transfer func¬ 
tion. In the PTA sensitive band, fcro ^ 0, then from the 

2(l-3it;) 

spectrum (8) we find oc A:”* 1 + 3 ™ . The obser¬ 

vations are influenced by the strength of RGWs, so they 
can merely determine the total tilt of ngw(A:), which is a 
combination of tit and w. 

Fig. 6 shows the spectrum of energy density undergoes 
the cosmic expansion and different EoS. For the unusual 
EoS with li; > 1/3, the energy density spectrum of RGWs 
is strongly intensified, while the matter-like EoS {w = 0) 
damps the RGWs greatly. It is clear that, when general 
EoS is considered, the spectrum at higher frequency is 
amplified or damped more than that at lower frequency. 

Tab. I lists the exact strength of flgw at three different 
frequencies. For the extremely stiff EoS with w = 00 
{w = 1) , the strength could be amplified at least 10^ 
(10^) times for the sensitive band of PTAs. For the mild 
EoS w = 0.6 [12, 59], Ggw could be amplified 9 times at 
10“® Hz, and as far as 130 times at 10“^ Hz. While a 
matter-like EoS {w = 0) may damp the strength by four 
orders of magnitude. 


w 

Hgw(lO-’^Hz) 

Hgw(10-®Hz) 

Hgw(10 ^Hz) 

00 

2.6 X 10"^® 

2.5 X 10"“ 

2.5 X 10“® 

1 

5.0 X 10"^® 

5.0 X 10”^'^ 

5.0 X 10"^® 

0.6 

9.2 X 10"^® 

3.4 X 10"^® 

1.3 X 10"^"^ 

1/3 

9.8 X 10"^'^ 

9.8 X 10”^'^ 

9.8 X 10"^'^ 

0 

3.8 X 10"^° 

3.8 X 10"^^ 

3.8 X 10"^"^ 


where j3 is related to the EoS parameter w through 

P = with w = -. (16) 

1 -I- 3w p 

For the era of radiation domination, w = ^ followed by 
a oc T, and for the era of matter domination, w = 0 
and a (X T^. Since the EoS before BBN (~ IMeV) is 
still not clear, in the general scenario, we set w as a free 
parameter, which represents the effective average EoS in 
this era. 

The transfer function caused by general EoS is [26] 


/ flb ' 


Va(ro)/ 

' \ k J 


(17) 


where Ob is the scale factor at BBN epoch (here we chose 
the epoch at T = IMeV without losing generality [2]) and 
H\y is the Hubble parameter at BBN epoch. The lower 
boundary fc > fcb is determined by the fact that only 
those modes entered before Tb are affected by unusual 
EoS. 

Note that the sensitive band of PTA, i.e. 10“® ~ 10“^ 
Hz, corresponds to 52 ~ 3.9 x 10^ MeV, which is above 
the lower boundary of BBN energy scale, therefore un¬ 
usual EoS could affect the detection of PTAs. In this 
sensitive frequency range, the tensorial index rij and the 
EoS w are degenerate. To see this clearly, we put the 


TABLE I: The strength of Hgw(/, tq) at / = 10“®, 10“® and 
10“^ Hz for EoS with w = 00 , 1,0.6, 1/3 and 0 respectively. In 
the calculation, we use r = 0.1 and nt = 0. The cosmological 
parameters are the same as those mentioned at the beginning 
of this section. 


III. CONSTRAINTS ON RGWS BY CURRENT 
PULSAR TIMING ARRAYS 

From Fig. 5 and Fig. 6 , we have seen that many damp¬ 
ing effects such as cosmological redshift, phase transition 
of QCD, the free-streaming dark fluids and the unusual 
EoS are well in the sensitive band of PTAs. Thus PTAs 
provide the good probe into these interesting physics. 
Now let us turn to the detection of RGWs through PTAs. 

RGWs could weakly stretch or suppress the space, 
leading to the change of time of arrival for electro¬ 
magnetic pulses from pulsars. Difference of time of ar¬ 
rival could be accumulated to form considerable timing 
residuals. By analyzing the timing residuals carefully it 
is possible to detect RGWs. 

In practice, the characteristic strain hc{f) is more con¬ 
venient for analysis and can be related to Dgw by [60] 

f^gw(/,ro) = 


(18) 
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FIG. 6: The dependence of f2gw(/, tq) on the EoS before BBN. 
Five kinds of EoS, including w = oo, 1, 0.6, 1/3 and 0, are 
considered. Here we choose r = 0.1, nt = 0 and adopt the 
cosmological parameters mentioned before. As the effect of 
unusual EoS changes the spectrum by several orders of magni¬ 
tude, we do not include the relatively small effects caused by 
cosmic phase transitions TpTik,To) and free-streaming parti¬ 
cles Tv{k, To) here, but the damping effect of cosmic expansion 
Tz(k,To) is still incorporated into the total transfer function. 
Sensitive band of PTAs is the shaded area in the figure. 


So if the spectrum of energy density (8) is given, the 
characteristic strain can be determined theoretically. 

When the spectrum is unknown, it is helpful to postu¬ 
late a parameterized power-law form as follows [23, 24], 

KU)=a{{X. (19) 


where A is the amplitude of hc{f) at /yr = 1/yr = 3.171 x 
10“® Hz, and a is the power index. 

Fig. 7 shows the current constraints on A from PPTA 
[22], EPTA [23] and NANOGrav [24]. From Fig. 7, we see 
that the result of PPTA is the most stringent one (about 
a factor of 1.5 better than NANOGrav at a = —2/3), but 
there is only one data point (the black diamond in the 
figure), so in the rest part of this paper we would con¬ 
centrate on results from NANOGrav and EPTA to cover 
more possible indices. In a wide range of a, NANOGrav 
gives a tighter constraint than EPTA for all a G [—2, 0] 
and the advantage becomes obvious as a increases. For 
example, at a = —1, NANOGrav gives A < 8.1 x 10“^®, 
which is only a factor of 1.7 better than 1.4 x 10“^^ 
of EPTA. When a increases to —0.6, NANOGrav gives 
A < 1.6 X 10“^®, a factor of 2.2 better than 3.5 x 10“^^ 
of EPTA. 


FIG. 7: The current 2(7-constraints on strain amplitude A 
from three PTA projects: NANOGrav, EPTA and PPTA. 
The NANOGrav data are obtained through figure 3 in [24]. 
The two red squares (a = —1,A = 1.4 x 10~^®) and {a = 
— 2f3,A = 3 X 10“^®) are constraints from EPTA by use of 
Bayesian method [23]. The dotted line is obtained by linearly 
fitting the squares. The black diamond (a = —2/3, A = 1.0 x 
10"^®) comes from PPTA [22]. 


A. Constraints on RGWs by current PTAs 
considering cosmic expansion only 

When upper limits on A for different indices are given, 
we can convert hdf) into Ogw(/) through (18) and com¬ 
pare it with the theoretical results in (8) to give limits 
on the inflationary parameters r and n*. 

We should firstly find out the relation between the two 
power indices a and n*. Because different damping ef¬ 
fects alter the spectrum differently (see Fig. 5), the rela¬ 
tion depends on the damping effects. Here we consider 
the spectrum damped by cosmic expansion only. In this 
case, the transfer function (7) has only one term T/(fc, tq) 
from cosmic expansion (10). For waves with frequencies 
27r/ ^ fceq, we find the simple relation between a and nt 

nt = 2a + 2, (20) 

and the relation between A and r 

327r^A^fcgqTo^ /fcpN^Y 1 
45Afl;(fco)f^m ^27r/ \fyj 

Inserting the parameters in Section H and the data in 
Fig. 7, we obtain the upper limits on r for different nt- 

Fig. 8 shows the constraints on r — nt space. The red 
dashed and the blue dashed line are the cases consider¬ 
ing Tz{k,To) only. At nt = 0, the limit from NANOGrav 
is 8.5 X 10®, which is a factor of 3.1 more stringent than 
2.7 X 10® from EPTA. Note that [23] gives a limit 2.5 x 10® 
for the same EPTA data. The small deviation occurs be¬ 
cause we adopt new cosmological parameters, which are 
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FIG. 8: 2(T-constraints on the r — nt space by NANOGrav [24] 
and EPTA[23]. The solid lines are the constraints for the case 
including transfer function of cosmic expansion Tz{k, tq) and 
that of cosmic phase transitions TpT{k,To), while the dotted 
lines contain Tz{k,To) only. 

slightly different from those used in [23]. As rit increases, 
the constraints become more stringent. At rit = 0.8, 
the constraint of r from NANOGrav is 0.04, while upper 
limit of EPTA is 0.17 (see Table II). For rit < 0 both 
constraints become very loose. In Table II, we also list 
the results for rit = —0.4. 

For r = 0.1, when only redshift is considered, Zhao et 
al [25] found rit = 0.90 for NANOGrav [61] and rit = 0.88 
for FPTA [62], based on the previous observations. While 
we found that, for r = 0.1, latest NANOGrav data [24] 
follows rit = 0.76 and EPTA data [23] follows nt = 0.83. 
Both constraints become tighter due to the update of 
NANOGrav and EPTA data. 


B. Constraints on RGWs by current PTAs 
considering phase transitions and neutrinos 

Although the damping effects beyond cosmic expan¬ 
sion Tz{f, To), such as cosmic phase transitions TpT'{k, tq) 
and the relativistic free-streaming gases Ti^(fc,To), have 
been well studied before, few authors include these damp¬ 
ing effects when they constrain the inflationary param¬ 
eters (see [25] for example). As we mentioned be¬ 
fore, cosmic phase transitions could damp the spectrum 
Ogw(/, Tq) as much as 60% in the sensitive band of PTAs 
(see Fig. 5). Thus, to make accurate detection, it is in¬ 
evitable to consider these damping effects. In this sub¬ 
section, we will discuss the constraints on the inflation¬ 
ary parameters r and rit when additional damping effects 
TpT(fc, To) and T^{k, tq) are included in the transfer func¬ 
tion. 

There are two points that deserve our attention. 
Firstly, it is very clear that neutrinos could not affect any 


detection approached by PTAs, as neutrinos only damp 
the spectrum in the range of 10“^® ^ 10“^® Hz, which 
is well beyond the detecting ability of PTAs. Secondly, 
TpT(k,TQ) makes the spectrum deviate from power law, 
so, strictly speaking, the postulation (19) does not hold. 
But TpT(k,To) is on the order of 0.1, which could only 
change the spectrum by a small amount, so in the mean¬ 
ing of perturbation, we could still apply the power law 
and attribute the small deviation from power law to the 
weak dependence of A on the frequency. Then the rela¬ 
tion (20) between a and tit remains, while the expression 
of r should be modified. 

When we include cosmic phase transitions, the transfer 
function T(fc,ro) will have another term TpT(fc,To) (11) 
besides Tz{k,To) (10), so the expression of r becomes 

^ ^ (221 

45A_R(fco)H^ V27r/ V/yr^ T|.p(fc,ro) 

As TpT(fc, To) < 1, we can expect the upper limit on 
r for fixed nt would increase after we include cosmic 
phase transitions. Specifically, numerical calculations 
give Tprp(fcyi., To) = 0.428, then we could obtain a looser 
constraint on r — nt space and the upper limit on r in¬ 
creases a factor of 2.3 for all indices. 

Fig. 8 shows the upper limit on r — nt- The red solid 
and the blue solid line are upper limits where Tp^{k, tq) 
is included. At nt = 0, NANOGrav [24] gives upper limit 
2.0 X 10®, and EPTA [23] gives 6.4 x 10®. Both are indeed 
raised by a factor of 2.3 compared to the limit without 
TpT(fc,To). At nt = 0.8, the constraint from NANOGrav 
is r < 0.09 and that given by EPTA is r < 0.39. We also 
list the results at nt = —0.4 in Table II. 


nt 

Include Tz only 

Include Tz and Tpt 

T nano 

'f' epta 

T nano 

'f' epta 

-0.4 

4.1 X 10® 

1.1 X 10^® 

9.6 X 10® 

2.6 X 10^® 

0 

8.5 X 10® 

2.7 X 10® 

2.0 X 10® 

6.4 X 10® 

0.8 

0.04 

0.17 

0.09 

0.39 


TABLE II: Constraints on upper limit of r for different nt by 
use of NANOGrav [24] and EPTA [23]. The second and third 
column are the limits when cosmic expansion Tz{k,To) is the 
only transfer function, while in the last two columns cosmic 
phase transition TpT{k,To) is added. 


C. Constraints on RGWs from other observations 

1. Constraints from BBN 

Besides PTA projects, BBN can also give upper limit 
on r — nt space. BBN precisely predicts the abundance 
of light elements in the Universe. Nucleosynthesis hap¬ 
pens efficiently when the collision rate of protons and 
neutrons are much higher than the Hubble parameter at 
that moment. To obtain enough abundance of H and 
He, the expansion rate of the Universe should not be too 
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large, otherwise the plasma soup would be diluted too 
much to produce effective reaction. Expansion rate of the 
Universe or Hubble parameter H depends on the energy 
density of radiation at that moment, thus the abundance 
of light elements could provide constraints on the total 
energy of radiation including RGWs. 

By use of the effective number of neutrino species IVeff, 
BBN could set an upper limit on the energy of RGWs [60] 



Hgw(fc,To)h^ d\nk < 5.6 x 10 ^{Ngg 


N,), (23) 


where Ni, is the number of neutrino species, while A^eff 
is the effective number of neutrino species. For stan¬ 
dard model in particle physics, considering the non- 
instantaneous decoupling of neutrinos, N^, = 3.046 [55]. 
Recent combined observations from BBN and primor¬ 
dial mass fraction of ^He and the abundance of deu¬ 
terium give iVeff = 3.41 with 95% confidence [57], then 
f dlnfc < 2.04 X 10“®. The lower boundary of 

the integration is determined by counting the RGWs 
that entered horizon at the BBN epoch only, while the 
upper boundary is determined by quantum limit. For 
the very short wavelength (i.e. high frequency) portion, 
the ultraviolet divergences is avoided by considering the 
Parker’s adiabatic theorem [63], which states that, dur¬ 
ing a transition between expansion epochs with a char¬ 
acteristic time during At, the gravitons created will be 
suppressed for wavenumbers k > 1/At. Here we follow 
[64] to choose /up = kup/^n = 10^° Hz by assuming the 
energy scale for the inflation is around 10^® GeV and 
fciow = k{T = IMeV) which can be obtained from (12). 

Inserting the transfer function T^ik, tq) (10) caused by 
cosmic expansion and TpT(fc, tq) by cosmic phase transi¬ 
tions (11) into the spectrum of energy density l}g^{k, tq) 
in (8), we can integrate the spectrum numerically to ob¬ 
tain the upper limit on the inflationary parameter space 
r — rit- The result is shown by the purple dashed line in 
Fig. 9. 


2. Constraints from CMB and matter power spectrum 

Due to the precise observations of temperature and 
polarization anisotropies, in the modern cosmology, 
GMB becomes one of the most important tools to con¬ 
strain various cosmological phenomena/processes includ¬ 
ing RGWs. RGWs affect the CMB fluctuations by two 
different ways. 

Firstly, RGWs, as well as density perturbations, are 
the sources which generate the CMB anisotropy power 
spectra, including the temperature anisotropy auto¬ 
correlation spectrum (TT), E-mode polarization auto¬ 
correlation spectrum (EE), B-mode polarization auto¬ 
correlation spectrum (BB), and the cross-correlation of 
temperature anisotropy and E-mode polarization (TE). 
Although the definite signal of RGWs has not yet been 
found in the CMB power spectra, based on the recent 



Ht 


FIG. 9: Constraints on the inflationary parameter space r — nt 
by use of data from observations of CMB (Planck [7] and the 
joint analysis of CMB and LSS (matter power spectrum) [56]), 
BBN [57], LIGO & Virgo [12], EPTA and NANOGrav [23, 24]. 
Note that the line stands for BBN is only a bit higher than 
that for GMB & LSS, so we use purple dashed line to represent 
the result of BBN and orange solid line for CMB & LSS. Here 
all the constraints have considered Tz{k,To), TpT{k,To) and 
T^{k, To). 


observations of TT, TE and EE data, Planck team gives 
the quite tight constraint on the tensor-to-scalar ratio 
r < 0.11 [7], which is nearly independent of the spec¬ 
tral index n* [65]. A similar constraint r < 0.12 is also 
obtained from the joint analysis of BICEP2/Keck Array 
and Planck B-mode polarization data [10]. 

Secondly, RGWs affect the growth of density pertur¬ 
bation and the CMB anisotropy power spectra by con¬ 
tributing extra energy density to the total energy density 
or, equivalently, by changing the expansion rate of the 
Universe. The weak interaction of RGWs with matter 
makes them similar to massless neutrinos. So constraints 
on additional neutrinos from CMB and large-scale struc¬ 
ture (LSS) observations could be applied to RGWs [66]. 
With 95% confidence level, CMB combined with matter 
power spectra gives Ngg — < 0.33 [56]. Again, we can 

insert this number into (23) to yield another constraints 
on RGWs. We obtain J dlnfc < 1.85 x 10“®, but 

the lower boundary of integration is much lower than 
that of BBN. We adopt /low = 10“^^ Hz [66] to do the 
integration. The corresponding constraint on r — rit is 
shown in Fig. 9 (the orange solid line). 


3. Constraints from gravitational-wave detectors 

The ground-based interferometer could directly detect 
RGWs in the range 10^ ~ 10^ Hz by analyzing the inter¬ 
ference stripes caused by gravitational waves in the laser 
beams. The latest constraints on the stochastic back¬ 
ground is from joint observations of LIGO and Virgo. In 
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frequency range 41.5 — 169.25 Hz, LIGO and Virgo cor¬ 
relations give an upper limit on energy density spectrum 
Hgw(fc,ro) < 5.6 X 10“® for rit = 0, while in the range 
600— 1000 Hz r2gw(fc,ro) < 0.14 for rit = 3 (there are the 
other two sets of constraints at nt = 0 and 3, but they 
are much looser, so we ignore them). Both limits are of 
95% confidence [12]. 

We can use the LIGO & Virgo data to estimate the 
upper limits of inflationary parameters r and rit- Firstly, 
we convert the limits on Hgw(A:,To) at a = 0 and 3 into 
strain amplitude A through (18), then use (22) to obtain 
limits on r and nt at these two indices. Finally, we lin¬ 
early fit the two points to cover the range in nt € (0,1), 
which is shown in Fig. 9 (the black solid line). 


4- Comparison with different constraints 

Fig. 9 shows the constraints on r — n* space from 
different observations. We can see that at low index 
(nt < 0.37), upper limit of r from Planck is the most 
stringent and constant constraint (r < 0.11), while at 
high index {nt > 0.37), constraint by use of joint anal¬ 
ysis of GMB and LSS becomes the most stringent and 
the upper limit decreases as nt increases, for example, at 
nt = 0.5, the upper limit is 10“^, but at nt = 0.8, the 
limit is r < 10“^^. 

We see that the constraints from BBN (the purple dot¬ 
ted line) is slightly looser than that from the joint data 
of GMB and LSS (the orange solid line). There are two 
reasons for it. Firstly, constraints on Ngg from BBN is a 
bit looser than that from the joint analysis (0.36 versus 
0.33, see subsection IB C 1 and IB G 2.) Secondly, the in¬ 
tegration range of BBN (/low 10“^° Hz) is larger than 
that of the combination of GMB and LSS (/low ~ 10“^^ 
Hz), thus less energy could be added to the integration, 
making the tensor index larger. 

Although the direct constraints from EPTA and 
NANOGrav have significantly improved recently, they 
are still much weaker than those from GMB, matter 
power spectra and BBN. So are LIGO and Virgo. There¬ 
fore, the great advancement is still needed to put tighter 
constraints. Advanced LIGO and future PTAs like FAST 
and SKA are indeed quite necessary. 


IV. DETECTION OF ROWS BY THE FUTURE 
OBSERVATIONS 

In this section we will forecast the detections of RGWs 
by the potential PTA observations, including the FAST 
and SKA projects. In addition, as we know, the damping 
effects caused by cosmological redshift, phase transition 
of QGD, the free-streaming dark fluid and the unusual 
EoS are all well in the sensitive band of PTAs. So, the 
future PTAs also provide a good chance to study these 
interesting physics, which will also be investigated in this 
section. 


In order to estimate detection abilities of the future 
PTAs, we will use time of arrival to establish the signal- 
to-noise ratio (SNR) of observations. When the spectrum 
of energy density Bgw (8) is given, an important prob¬ 
lem follows: What confidence level can we achieve when 
we come up with new PTA projects to detect RGWs? 
The key lies on the timing residuals caused by RGWs 
and their correlations. RGWs could stretch or suppress 
the space weakly, leading to the change of time of arrival 
for electro-magnetic pulses from pulsars. Difference of 
time of arrival could be accumulated to form consider¬ 
able timing residuals. The timing residuals are usually 
decomposed into two parts 

R^{tk) = s\tk) + n\tk), (24) 


where R''{tk) is the timing residual for i-th pulsar at time 
tfe, while s'‘{tk) and n’’{tk) are the timing residuals con¬ 
tributed by RGWs and by noise. To make a detection, it 
is necessary to discriminate the contribution of RGW sig¬ 
nal from that of noise. Let us consider the characteristic 
of s®(ffc) and n'‘{tk) separately. 

The genuine RGWs would induce common residuals in 
the data, so the correlation of s®(ffc) gives [67] 

{s\tk)s^{tk'))=alH{e,,)5t,,t,,, (25) 

where St^^ty is Kronecker delta function and H{6ij) is 
the Hellings-Downs curve, while cTg is root mean square 
(RMS) of timing residuals resulted from RGWs. The 
Hellings-Downs curve [20, 68] is 

H{0^j) = '^x\nx-^ + ]^{l+5{x)^,x= - —(26) 


where 9ij is the angular separation between z-th and j- 
th pulsar. The RMS of timing residuals from RGWs is 
given by 


(7 


2 

g 


y'‘ \hc{f)\ 

Jf, 


df, 


(27) 


where // = 1 /T and fh = 1/{2At) are the lower and the 
higher boundary of sensitive band. T and At are the 
span of whole observation and the interval between two 
observations respectively. Now, if the spectrum of energy 
density (8) is given, the RMS of RGW is determined. 

As for the part of noise, we assume that they are white 
and the same for every pulsar in the PTA, then the cor¬ 
relation of contributions from the noise is simple [67] 

{n\tk)n^ {tk')) = crlSijStyty, (28) 


where is the RMS of timing residuals from white noise. 

By use of the correlation of timing residuals and the 
whitening method, we could obtain the SNR for PTA 
[25, 69, 70] 


(SNR) = Vn\ 1 -h 


E. 


1-b 


VTVa 




Pd(A')J 


(29) 
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where brackets ( ) are ensemble average (hereafter when 
we say SNR we mean the expectation of SNR, or (SNR)). 
N = n{n — l)/2 is the number of pairs with n be¬ 
ing the number of pulsars. H'^ and are the aver¬ 
age of H‘^{9ij) and the variance of H{6ij) respectively. 
For pulsars evenly distributed in the full sky, we have 
W=Y'\j = 1/48. 

As for Pg(A'), it is defined as the RMS in the A'-th 
bin [69]: 


A'+0-5 


^g(A') = 


A^-0.5 


hlU) 

1271^/3 


d/, (A' > 1) (30) 


where A' is the bin A whose signal is higher than noise, 
i.e. Rg(A) > cr^/m, with m being the total number of 
observations in the PTA project. For A' = 1 the lower 
boundary of integration in Pg(A') is 0.97/T [69]. The 
total spectral density in a bin is the sum of contributions 
from signals and noises: Pd(A') = Rg(A') + a^jm. 

We could see that, to calculate the SNR for a given 
spectrum of energy density, an observational plan should 
at least include the following information: the observa¬ 
tion span T, the number of pulsars n, the number of 
observations m and the noise level. 

For timing noise, due to the intrinsic instability in the 
pulsars and uncertainties in the telescopes, the noise is 
not guaranteed to be white and may have different value 
for different pulsars [67]. In our analysis, for simplicity, 
we assumed white and equal noise level Un for all pulsars 
in a given array. 


A. Potential PTA projects in the future 

In this paper, we follow Zhao et al [25] to consider 
the following future observations: the complete PPTA, 
FAST, SKA and optimal PTA. 

The complete PPTA plans to observe 20 pulsars which 
have the same noise level of 100 ns over 5 years (T = 5 
yr, tTn = 100 ns, n = 20, to = 250), while current PPTA 
mainly depends on four pulsars [22]. So the future PPTA 
could greatly improve the observation. 

Another observation plan that deserves our attention is 
that from FAST, which would be the largest single dish in 
the world when it obtains the first light in 2016. Working 
with multi beams in the 70 MHz—3GHz frequency band, 
FAST is expected to discover 4000 Galactic pulsars and 
one tenth may be millisecond pulsars [71]. For the search 
for gravitational waves, FAST plans to run a series of 
similar observations (T = 5 yr, n = 20, to = 250, we 
use FAST(20) to label this plan) but with much lower 
noise of only 30 ns [71], which may directly improve the 
value of SNR. Here we assume a second plan, FAST(40), 
to adequately cover the possibility of more pulsars in a 
longer duration of timing: T = 10 yr, n = 40, to = 500. 

In the 2020s, the establishment of the thousands of 
single dishes will erect the biggest radio array SKA and 


pulsar searching and timing will be part of the SKA sci¬ 
entific goals. Accounting for the large signal-collecting 
area, SKA could well survey the pulsars in the Milky 
Way in an unprecedented efficient way [72] and greatly 
facilitate the detection of gravitational waves. To esti¬ 
mate the capacity of SKA in detecting RGWs, we adopt a 
low-noise plan running for ten years (T = 10 yr, Un = 50 
ns, n = 100, TO = 500, we use SKA(IOO) to denote this 
plan) [73]. However, from a conservative view point, we 
would also consider a less ambitious one, SKA(40), with 
T = 10 yr, cTn = 50 ns, n = 40, to = 500. 

Prolonging time span allows PTAs to approach more 
low frequency signals, which can be much stronger than 
that of high frequency. So, accounting of these improve¬ 
ments, future PTAs may greatly increase their sensitiv¬ 
ities. Zhao et al [25] also considered the optimal PTA 
(T = 20 yr, a-^ = 30 ns, n = 200, to = 1000), which can 
reach an unprecedented sensitivity. For easy reference, 
we list these projects in Table HI. 


Potential PTA T/yr nn/ns n m 

Complete PPTA... 5 100 20 ^ 

FAST(20). 5 30 20 250 

FAST(40). 10 30 40 500 

SKA(40). 10 50 40 500 

SKA(IOO). 10 50 100 500 

Optimal PTA. 20 30 200 1000 


TABLE III: The parameters of potential PTA projects: Time 
span r, noise level a^, number of observed pulsars n and the 
total number of observations m. 


For the known damping factor, if r and rit are given, 
the spectrum (8) is determined. Then, we use (18) to find 
the characteristic strain hc{f) and obtain Pg(A') through 
(30) for known PTA observations. Finally, we can calcu¬ 
late the SNR for these observations through (29). 

Note that without RGW signals, the timing residuals 
(24) will be Gaussian distributed, thus SNR = 2 repre¬ 
sents 95% confidence level. In the calculation, we would 
like to obtain the constraints on r — nj space with 95% 
confidence for the potential observations. This is possi¬ 
ble, because for a given observational plans, when SNR 
is fixed, r and nt are constrained by the relation (29). 
So we firstly insert the spectrum with free parameters r 
and rit into (18) to obtain Pg(A') through (30), then set 
SNR = 2 to give upper limits on r for different rit with 
95% confidence level. 


B. Detecting RGWs in the accelerating Universe 

As the first step, we only consider transfer function 
due to cosmic expansion Tz{k,TQ) in (10) and assume 
a normal radiation EoS w = 1/3 before BBN in this 
section. 

Fig. 10 show the SNR of FAST in the range rit S [0,1]. 
The black and grey lines are the SNR including damp¬ 
ing effect from Tz{k,To) only. Our results for FAST(20) 
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FIG. 10: The SNR of FAST(20) and FAST(40) when differ¬ 
ent damping effects are considered. For FAST(40), the black 
lines consider the effect of cosmic expansion Tz{k,To) only, 
while the red ones include cosmic phase transitions TpT{k, to). 
For FAST(20), the grey lines consider Tz(k,To) only and the 
orange ones add TpT(k, tq). In each project, from up to 
down, the solid, dashed and dotted lines are for the case of 
r = 0.1, 0.01, 0.001 respectively. 


are the same as those in [25]. At rit = 0, SNR vanishes 
for both FAST plans, thus no RGWs can be detected. 
However, as nt increases, the SNR increases a lot. If we 
fix rit = 0.8, when r = 0.1, RGWs produce a signifi¬ 
cant detection with SNR= 4 for FAST(20) and 11.5 for 
FAST(40). 

It is obvious that smaller r leads to lower SNR. For 
example, at nt = 0.8, when r decreases from 0.1 to 0.001, 
SNR jumps from 11.5 to 6.5 for FAST(40), and from 4 
to 2 for FAST(20), thus greatly weaken the detection. 

The trends described above also hold for SKA, as we 
could see from Fig. 11, which shows the SNR for two 
SKA plans with 40 and 100 pulsars. Since the only dif¬ 
ference between SKA(40) and SKA(IOO) is the number 
of pulsars (see Table. Ill), we conclude that increasing 
the scale of arrays could greatly raise the detection abil¬ 
ity. Lower noise level also increases the possibility to 
make a detection, therefore, FAST(40) will performe bet¬ 
ter than SKA(40). For example, at nt = 0.8, for r = 0.1, 
FAST(40) gives SNR=11.5, which is a bit higher than 10 
given by SKA(40). 

In Fig. 12, we plot the constraints on r — nt space 
by potential observations. The dashed lines are the 
cases only considering Tz{k,To). The current limits from 
NANOGrav and EPTA are also plotted here to make 
comparison. We could see, when only Tz{k, tq) is consid¬ 
ered, in the whole range nt € [0,1], for the upper lim¬ 
its on r, complete PPTA could be about 10 times more 
stringent than current NANOGrav and about 50 times 
better than current EPTA, while EAST(20) is 11 times 
more stringent than complete PPTA and SKA(IOO) is 


FIG. 11: The SNR of SKA(40) and SKA(IOO) when T^{k,To) 
and TpT{k, tq) are considered. The color is similar to Fig. 10. 


137 times better than FAST(20). The most stringent 
limit would come from the optimal PTA, which is about 
107 times better than SKA(IOO). 

Although both SKA plans could defeat the FAST(20), 
the advantages of SKA would be greatly challenged by 
FAST(40), which is a factor of 2.5 better than SKA(40) 
and only a factor of 1.6 weaker than SKA(IOO). The 
good competence of FAST(40) comes from the low tim¬ 
ing noise, as indicated by Table. III. As SKA(IOO) and 
FAST(20) are the upper and lower limits of the four plans 
of SKA and FAST, we would concentrate on these two 
observations in the rest part of this paper. However, due 
to the similar resolution of FAST(40) and SKA(IOO) in 
the r — nt space, the results for SKA(IOO) also hold for 
FAST(40) with good accuracy. 

At n* = 0, upper limit on r given by optimal PTA is 
0.44, which is the best constraints. At nt = 0.8, upper 
limit given by complete PPTA is 0.008, while FAST(20) 
gives 0.0007 and SKA(IOO) and the optimal PTA are 
much more stringent. 

For r = 0.1, the optimal PTA gives nt < 0.08, while 
for SKA(IOO) the upper limit is 0.32, for FAST(20) it is 
0.58, and for complete PPTA, it is 0.69. If r = 0.01, the 
optimal gives limit nt < 0.18, SKA(IOO) gives nt < 0.44, 
and FAST(20) gives nt < 0.68. So both SKA(IOO) and 
optimal PTA could give fairy good constraints on the 
inflation models. 


C. Effects of phase transitions and neutrinos 

In the previous work [25] , we only calculated the SNR 
for the RGWs damped by cosmic expansion Tz{k,To), 
which is briefly reviewed in the previous subsection. 
However as we mentioned before, cosmic phase transi¬ 
tions, especially the phase transition of QGD, damp the 
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FIG. 12: The constraints on r — nt space with 95% con¬ 
fidence by future PTAs described in Table. Ill: Complete 
PPTA, FAST, SKA and Optimal PTA. Dashed lines are the 
constraints considering cosmic expansion Tz{k, to) only, while 
solid lines consider both Tz{k,Ta) and cosmic phase transi¬ 
tions TpT(fc,To). The lines of NANOGrav and EPTA are 
copied from Fig. 9 for comparison. Here, for clarity, we only 
plot the solid lines for FAST(40) and SKA(40). 


spectrum by amount of 60% at / > 10“® Hz, which is well 
in the sensitive band of PTAs. Relativistic free-streaming 
neutrinos could also damp the spectrum by about 40% in 
the range of 10“^® ~ 10“^® Hz. Although the behavior of 
neutrinos has no effect on PTAs, they play an important 
role in the constraints by use of combined CMB and mat¬ 
ter power spectra. Therefore effects from cosmic phase 
transitions and neutrinos deserve attention. 

In this subsection, we shall take into account the damp¬ 
ing effects caused by cosmic phase transitions TpT(fc,To) 
in (11) and neutrinos T^{k,To) in (13). Here we also as¬ 
sume a normal EoS with w = 1/3. 

For PTAs we only deal with the effects of TpT(fc,ro) 
due to QCD transition and ignore the damping effects of 
e“'"e“ annihilation and SUSY breaking, as e^e“ annihi¬ 
lation affects RGWs in the range of / < 1.7 x 10“^^ Hz 
(or, equivalently, 0.1 MeV) and SUSY phase transition 
in the range of / > 3.1 x 10“® Hz (or 10® MeV). We add 
the transfer function TpT{k, tq) and Ti,(k, tq) to the total 
transfer function, and follow the procedures described in 
Section IV A to obtain SNR and the constraints on r and 
rit respectively. 

In Fig. 10, the red and orange lines show the SNR 
for RGWs further damped by TpT(k,To) in FAST plans. 
Obviously, compared with the cases which only contain 
Tz{k,To), SNR decreases for both plans when additional 
cosmic phase transitions are considered. The decline 
of SNR is consistent with our intuition, as additional 
damping could surely lead to weaker RGW signals and 
smaller SNR. In each observation, when rit > 0.6, for 
r = 0.1, SNR decreases by an almost constant amount: 



rit 


FIG. 13: The constraints on r — nt space with 95% by future 
observations compared with current observations from CMB 
and BBN. Here, the PTAs consider both cosmic expansion 
Tz{k,To) and phase transitions rpT(fc,To). The lines of CMB 
and BBN are copied from Fig. 9, while those of PTAs from 
Fig. 12. 


FAST(40) decreases one unit, while FAST(20) by amount 
of 0.6. So cosmic phase transitions are more important 
in FAST(40) than in FAST(20). As nt decreases, the dif¬ 
ference in SNR caused by TpT(fc,ro) also decreases and 
tends to zero at n* = 0. Similar features appear in SKA 
observations too (see Fig. 11). 

Fig. 12 shows the effects of cosmic phase transitions 
TpT(fc, To) on parameters r and nt- Obviously, Tp'pik, tq) 
raises the upper limits of all PTAs nearly by a constant 
factor of 2.2 in the wide range of nt- It is not difficult to 
understand the reason, as the big jump cased by QCD 
transition is the main feature in the narrow sensitive band 
of PTA (see Fig. 5). 

In Fig. 13, we put the constraints from future PTAs 
and those from current observations together to get a 
better view. The bounds from PTAs have already in¬ 
cluded damping effects from cosmic expansion Tz{k,To) 
and phase transitions TpT(fc,To). We see that around 
nt = 0.38, SKA(IOO) could do much better in limiting r 
and become very competitive with current CMB observa¬ 
tions from Planck and the joint analysis from CMB and 
LSS. In nt € (0.36,0.38), SKA(IOO) is a bit more (about 
a factor of 1.3) stringent than BBN and CMB. In com¬ 
parison with constraints from CMB and BBN, the role of 
cosmic phase transitions becomes particularly important 
for SKA(IOO), because if without rpT(fc, tq) the advan¬ 
tage of SKA(IOO) could increase to about 3 times better 
than CMB and BBN over a wider range. We can also 
find that optimal PTA could replace the dominant role 
of CMB and BBN in the range n* € (0.11,0.5) and re¬ 
strict nt < 0.23 for r = 0.01 and nt < 0.34 for r = 0.001, 
which are very tight limits for a wide variety of inflation 
models, including slow-rolling inflation. We list the up- 
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per limits on rit in Table IV for clarity. When we fix 
rit, we see that at rit = 0, CMB is still most stringent 
constraints even for the optimal PTA, but for a small 
positive tensorial index, e.g. rit = 0.2, the optimal PTA 
can require r < 0.02. Therefore, the future constraints 
are much better than the present ones if the spectrum of 
RGWs is blue tilted (i.e. n* > 0). The bounds on r for 
different indices n* are listed in Table V. 


r 

upper limit 

on nt 

without TpT(k,To) 

with TpT{k,To) 

0.1 

0.08 

0.12 

0.01 

0.18 

0.23 

0.001 

0.30 

0.34 


TABLE IV: The best constraints on tensorial index nt with 
95% confidence level for three typical tensor-to-scalar ratios: 
r = 0.1,0.01 and 0.001. Here we choose a normal EoS with 
ui = 1/3. As the optimal PTA would be the best constraint, 
these bounds with cosmic phase transitions TpT{k,To) are all 
from optimal PTA in Fig. 13, while those without phase tran¬ 
sition are read out from Fig. 12. 



FIG. 14: The impact on SNR of SKA(IOO) caused by dark 
fluids decoupling at 10^, 10® and 10® MeV compared with the 
case without dark fluid. Here we choose r = 0.1 and pick out 
the part with obvious difference among these three cases. 


nt 

^ cmb&lss 

rska(lOO) 

^optimal pta 

0 

> 1 

> 1 

1 

0.2 

> 1 

> 1 

0.02 

0.4 

0.03 

0.04 

3.2 X 10“® 

0.6 

1.0 X 10“'^ 

6.3 X 10“® 

6.0 X 10“® 

0.8 

5.0 X 10“®® 

1.0 X 10“® 

1.0 X 10“'^ 


TABLE V: The constrains on the upper limit of r with 95% 
confidence for fixed nt from CMB & LSS, SKA(IOO) and the 
optimal PTA. Note that the upper limit given by pure CMB 
(Planck) observation is r < 0.11. All the constraints have con¬ 
sidered the damping effects from cosmic expansion Tz(fc, ro), 
cosmic phase transitions rpT(fc,'ro) and the free-streaming 
neutrinos Ti,{k,To). We use a normal EoS with w = 1/3 
here. 


D. Effects of relativistic free-streaming dark fluids 

The possible excess of effective number of neutrino 
species Nes leads to tension between observations and 
theories. As we mentioned in Sec. IIC 2, this problem 
may be reconciled by assuming extra relativistic light 
particles, such as massless bosons proposed by [58] or 
a class of massless fermions i.e. the dark fluid in this 
paper. If verified, the existence of extra species of parti¬ 
cles will bring new physics both to particle physics and 
cosmology. 

For RGWs, we have found that dark fluid plays dual 
roles in the evolution of the gravitational waves (see Fig. 
5): It damps the spectrum by adding the fraction of free- 
streaming particles (see Fig. 3). Meanwhile, it raises the 
RGW spectrum a bit through cosmic phase transitions 
(see Eq. (14)). 


In this subsection, we will focus on the impact of dark 
fluid on SNR and constraints on r — rit space. As before, 
we also assume an EoS with w = 1/3 before BBN. Since 
dark fluid could affect the spectrum through both cos¬ 
mic phase transition and free-streaming gases, we should 
combine the transfer function TpT'{k,To), T^{k,TQ) and 
Tz{k,TQ) to get the total transfer function in the spec¬ 
trum (8), then perform the similar procedures depicted 
in Sec. IV A to calculate SNR of RGWs for different po¬ 
tential observations. 

Fig. 14 shows the difference on SNR of SKA(IOO) 
caused by dark fluids decoupling at 10^, 10® and 10® 
MeV. It is not strange to find that dark fluid can raise 
the SNR curve a bit (< 0.3), as dark fluid could raise 
the spectrum through phase transition and counteract 
the damping effect caused by its free-streaming (see Eq. 
(14)). We could see this clearly in Eig. 5: For the dark 
fluid decoupling at 10 MeV (the black solid line, which 
is partially covered by the green solid line for 10^ MeV), 
the spectrum in the PTA band is clearly raised a bit, 
while for that decouples at 10® MeV (the blue dashed 
line), the spectrum is not only raised a bit at the high 
frequency end (i.e. / ~ 10“^ Hz), but also damped a bit 
at the frequencies below QGD transition (i.e. / ~ 10“® 
Hz). While at the middle frequencies (/ ^ 10“® Hz), 
the two effects nearly cancel each other. For dark fluids 
that decouple at 10® MeV, the effect of free-streaming 
could extend to higher frequencies, making the spectrum 
nearly unchanged. Therefore the SNR for the case of 10® 
and 10® MeV are raised more than that for 10® MeV (the 
spectrum for 10® MeV increases because the effect near 
10“®’ Hz become dominant as rit >0). This is consistent 
with Fig. 14, where the SNR curve for dark fluid decou¬ 
pling at 10® MeV (the purple solid line) is the middle one 
between the curve for 10® MeV (the green dotted line) 
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and that without dark fluid (the orange solid line). 

As rit goes up to 1, the spectrum could be raised more 
by dark fluids, because the contribution of high frequen¬ 
cies become dominant for n* > 0. We find the largest 
influence on SNR caused by dark fluid for n* > 0.9, but 
the SNR of SKA(IOO) is only raised by less than 0.5 for 
r = 0.1. Obviously, the impact on FAST would be even 
smaller (similar to the cases in Fig. 10). So the influence 
of dark fluid could reasonably be ignored by PTAs and 
we cannot expect any significant changes caused by dark 
fluid in the r — rit space. 

The reason why dark fluid plays small influence on the 
spectrum is simple. Current constraints on the number 
of species Ngg is so stringent that only one or two more 
kinds of extra particles are allowed, thus the contribution 
of dark fluid to the total radiation around the PTA band 
(10“® ~ 10~^Hz, or equivalently, 52 ^ 3.9 x 10^ MeV) is 
quite small, so both phase transition and free-streaming 
of dark fluid become insignificant. 


E. Effects of different equation of state w 

BBN could trace physics up to about 10 MeV in the 
history of Universe. However, the physics above this en¬ 
ergy scale is not quite clear, for example, EoS in this 
range is not determined. Reheating after inflation and 
the appearance of massive particles may lead to quite 
interesting EoS, which could make the Universe evolve 
differently. For RGWs, these unusual EoS may greatly 
amplify or depress the amplitudes, making the detection 
hopeful or despaired (see Fig. 6). With the advancement 
in the detection of RGWs, more stringent constraints 
could be set to the EoS in this early stage of the Uni¬ 
verse and limit the possible reheating physics. 

Since the transfer function of unusual EoS alters the 
RGW spectrum dramatically (up to several orders of 
magnitude), we will neglect the damping effects induced 
by cosmic phase transitions and free-streaming particles 
in this subsection, and only consider effects caused by 
cosmic expansion Tz{k,To) and unusual EoS T'Eos(fc) t'o)- 
Similar to the previous works [12, 26], we will consider 
five different EoS with w = oo, 1, 0.6, 1/3 and 0 to cover 
the stiff, soft and radiation-like states. 

Eirstly, we consider, when unusual EoS T'Eos(^i'd)) ap¬ 
pears, to what extent the current NANOGrav data can 
limit the r — rit space. With procedures similar to (21) 
and (22), we could find the expression of r. Then by use 
of the NANOGrav data (see Eig. 7), we can calculate the 
upper limits on r for the unusual EoS mentioned above. 
Eig. 15 shows the constraints on r — n* space when both 
Tz{k,To) and TEos(fc, tq) are considered. We can see that 
stiff EoS {w > 1/3) could greatly suppress the possible 
parameter space. Eor example when r = 0.1, the upper 
limit on rit is 0.6 for w = 0.6 and 0.49 for w = 1, both of 
which are more compact than that of 0.76 for w = 1/3. 
The upper limit could even be pushed to 0.21 for w = oo. 

Now we consider the potential constraints from the 



FIG. 15: The constraints on r and nt with 95% confidence 
level by use of NANOGrav data [24]. Here the transfer func¬ 
tion of cosmic expansion Tz{k,TQ) and that of unusual EoS 
TEosik,To) are considered. Note that the constraints of EoS 
with w = 0 are quite loose and exceed the scale of this plot, 
while the line for ui = oo is linearly extended to nt = 1. 


future PTAs. Inserting the transfer function TEos{k,To) 
of unusual EoS and Tz{k,To) of cosmic expansion to the 
total transfer function in the RGW spectrum, we could 
follow the steps given in Sec. IV A to calculate the SNR 
and the constraint on r — n* for different observations. 



FIG. 16: The SNR of SKA(IOO) when different EoS before 
BBN are considered. Here we choose r = 0.1 and plot SNR 
for five different EoS. The damping effect of cosmic expansion 
Tz{k,To) is also included in the calculations. 

Fig. 16 shows the SNR for different EoS in the 
SKA(IOO) project. We can see that, compared with the 
normal EoS of radiation {w = 1/3), stiff EoS {w > 1/3) 
could greatly increase the SNR of PTAs in the whole 
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range of rit £ (0.15,1). SNR even starts to saturate 
with SNR = 67.1 at rit = 0.6 for the inhnity stiff EoS, 
and reaches maximum at nt = 1 for w = 1. Around 
rif = 0, SNR for w = 1 and infinity can reach a well 
level. For example, when rit = 0.07, we have SNR=2 
for w = 1. Therefore, SKA(IOO) is still very sensitive to 
the EoS above li; > 1. If the EoS before BBN is dom¬ 
inated by matter of w = 0, then the SNR will be dra¬ 
matically damped in the whole range of n* £ (0,1) and 
the maximum of SNR is 8.9 at n* = 1. So, in conclusion, 
RGWs are quite sensitive to the EoS before BBN and 
could provide good probes into the interesting physics in 
this epoch. 

Fig. 17 shows the constraints on r — n* for SKA(IOO) 
when different EoSs are considered. We could see that, 
stiff EoS with w > 1/3 makes the constraints more 
stringent than the case of re = 1/3, while the matter-like 
EoS {w = 0) damps the spectrum so much that the 
constraint becomes very loose. Explicitly, at r = 0.1, the 
upper limit on nt given by w = oo is less than 0, while 
w = 1 gives 0.08, and w = 0.6 gives 0.18, all of which are 
more stringent than 0.32 given by normal radiation-like 
EoS {w = 1/3), but if w = 0, the upper limit on nt is 
near 1. So we can conclude that stiff EoS could facilitate 
the constraints of RGW parameters. For the case of 
r = 0.1, 0.01 and 0.001, we list the constraints on nt in 
Table VI. 


w 


upper limit of nt 


r = 0.1 

r = 0.01 

r = 0.001 

oo 

< 0 

< 0 

0.04 

1 

0.08 

0.18 

0.29 

0.6 

0.18 

0.29 

0.40 

1/3 

0.32 

0.44 

0.55 

0 

0.84 

0.94 

> 1 

TABLE VI: 

Upper limit 

of nt for different 

EoS given by 


SKA(IOO) with 95% confidence level when r is chosen to be 
the three typical values: r = 0.1, 0.01 and 0.001. The damping 
effect from cosmic expansion Tz{k,To) is also included here. 


V. CONCLUSIONS 

A stochastic background of relic gravitational waves, 
generated during the early inflationary stage, is a ne¬ 
cessity dictated by general relativity and quantum me¬ 
chanics. The spectrum of RGWs directly depends on 
the inflationary physics. So, it is always treated as the 
smoking-gun evidence of inflation. In addition, in the 
post-inflation stage, various cosmic phase transitions and 
relativistic free-steaming gases, also left imprints on the 
RGW spectrum. For this reason, RGWs also provide the 
cleanest way to probe the physics in the post-inflation 
epoch. 

Detections of RGWs have been carried on by different 
experimental methods. In the median frequency range 



FIG. 17: Constraints on r — nt space with 95% confidence 
level in SKA(IOO) for different EoS {w = oo, 1,0.6,1/3 and 
0) before BBN epoch. Here we consider damping effects 
from cosmic expansion Tz{k,To) and that from unusual EoS 
TEos{k,To) before BBN. 

/ £ (10“®, 10“^) Hz, the detection is by analyzing the 
timing residual of the millisecond pulsars. Recently, all 
three PTA groups (PPTA, EPTA and NANOGrav) re¬ 
ported the latest constraints on the GW background. In 
this paper, by considering the current and the potential 
future PTA observations, we investigated the constraints 
on the RGWs and the inflationary parameters in the gen¬ 
eral cosmological scenario. In particular, we studied the 
effects of cosmic phase transitions and various relativistic 
free-streaming fluids, which had been neglected in all the 
previous works. 

Gosmic phase transitions, including e“'“e“ annihilation, 
QGD transition and SUSY breaking, damp the RGW 
spectrum in the frequency range / > 10“^*^Hz, which is 
exactly the sensitivity range of PTA method. Taking into 
account the damping effects caused by all physical transi¬ 
tions, we find the upper limit of the tensor-to-scalar ratio 
r increases by a factor ~ 2 for any given spectral index 
Uf. In the standard cosmological scenario, for current 
NANOGrav constraints with nt = 0, the upper limit of r 
increases from 8.5 x 10® to 2.0 x 10®. While, for the future 
SKA(IOO), if r = 0.1, we find the detection of RGWs is 
possible only if nt > 0.36, instead of nt > 0.32. 

The relativistic free-streaming gases, including neutri¬ 
nos and some unknown dark fluids in the Universe, influ¬ 
enced the evolution of RGWs in the radiation-dominant 
stage, which significantly damped the RGWs spectrum at 
the frequency range / > 10“^®Hz. By analysis, we find 
this effect is very small in frequency / £ (10“®, 10“^) 
Hz for both neutrinos and dark fluids. Even for the fu¬ 
ture SKA(IOO) project, this effect seems impossible to be 
detected. 

In the general cosmological scenario, the effective EoS 
w of the cosmic fluid can be different from 1/3 in the 
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pre-BBN stage. We find that, the value of w greatly 
affects the RGWs spectrum, as well as their detection by 
pulsar timing arrays. For the future SKA(IOO) project, 
if r = 0.01 is set, we find the detection is possible if rit > 
0.44 for the standard model with w = 1/3. However, if 
w = 1, they can be detected only if m > 0.18. So, the 
stiff EoS could significantly decreases the difficulties of 
RGW detection. 

At the end of this paper, we should mention that 
in addition to the detection methods above, RGWs 
have also been constrained by some other observational 
and experimental efforts. In the frequency range / € 
(10“^®, 10“®)Hz, GW background produces a pattern of 
apparent proper motion of quasars. So, by observing the 
motion of quasars in the Universe, an interesting con¬ 
straint Ugw ^ 0.2 in this frequency band is given in [74]. 
For the RGWs with high frequency / > 10® Hz, by ana¬ 
lyzing the implications of graviton to photon conversion 


in the presence of large-scale magnetic fields, an upper 
limit Hgw ^ 1 is derived in [75]. Meanwhile, the experi¬ 
mental detections of the RGWs by various GW detectors 
(e.g. the cryogenic resonant bar detectors [76], the cav¬ 
ity detectors MAGO[77], the waveguide detectors [78], 
the Gaussian maser beam detectors [79]) have also been 
well studied in the recent literatures. 
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